{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "import matplotlib.pyplot as plt "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    x1, x2 = x\n",
    "    return 2 * x1 ** 2 - 2 * x1 * x2 + x2 ** 2 + 2 * x1 - 2 * x2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def grad(f,x): \n",
    "    h = np.cbrt(np.finfo(float).eps)\n",
    "    d = len(x)\n",
    "    nabla = np.zeros(d)\n",
    "    for i in range(d): \n",
    "        x_for = np.copy(x) \n",
    "        x_back = np.copy(x)\n",
    "        x_for[i] += h \n",
    "        x_back[i] -= h \n",
    "        nabla[i] = (f(x_for) - f(x_back))/(2*h) \n",
    "    return nabla \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.055454452393339e-06"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.cbrt(np.finfo(float).eps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def line_search(f,x,p,nabla):\n",
    "    '''\n",
    "    BACKTRACK LINE SEARCH WITH WOLFE CONDITIONS\n",
    "    '''\n",
    "    a = 1\n",
    "    c1 = 1e-4 \n",
    "    c2 = 0.9 \n",
    "    fx = f(x)\n",
    "    x_new = x + a * p \n",
    "    nabla_new = grad(f,x_new)\n",
    "    while f(x_new) >= fx + (c1*a*nabla.T@p) or nabla_new.T@p <= c2*nabla.T@p : \n",
    "        a *= 0.5\n",
    "        x_new = x + a * p \n",
    "        nabla_new = grad(f,x_new)\n",
    "    return a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "def BFGS(f,x0,max_it):\n",
    "    '''\n",
    "    DESCRIPTION\n",
    "    BFGS Quasi-Newton Method, implemented as described in Nocedal:\n",
    "    Numerical Optimisation.\n",
    "    INPUTS:\n",
    "    f:      function to be optimised \n",
    "    x0:     intial guess\n",
    "    max_it: maximum iterations \n",
    "    plot:   if the problem is 2 dimensional, returns \n",
    "            a trajectory plot of the optimisation scheme.\n",
    "    OUTPUTS: \n",
    "    x:      the optimal solution of the function f \n",
    "    '''\n",
    "    d = len(x0) # dimension of problem \n",
    "    nabla = grad(f,x0) # initial gradient \n",
    "    print(nabla)\n",
    "    H = np.eye(d) # initial hessian\n",
    "    x = x0[:]\n",
    "    it = 2 \n",
    "    while np.linalg.norm(nabla) > 1e-5: # while gradient is positive\n",
    "        if it > max_it: \n",
    "            print('Maximum iterations reached!')\n",
    "            break\n",
    "        p = -H@nabla # search direction (Newton Method)\n",
    "        print(p)\n",
    "        # p = np.linalg.solve(H, nabla)\n",
    "        a = line_search(f,x,p,nabla) # line search \n",
    "        print(a)\n",
    "        s = a * p \n",
    "        x_new = x + a * p \n",
    "        nabla_new = grad(f,x_new)\n",
    "        y = nabla_new - nabla \n",
    "        y = np.array([y])\n",
    "        s = np.array([s])\n",
    "        y = np.reshape(y,(d,1))\n",
    "        s = np.reshape(s,(d,1))\n",
    "        r = 1/(y.T@s)\n",
    "        li = (np.eye(d)-(r*((s@(y.T)))))\n",
    "        ri = (np.eye(d)-(r*((y@(s.T)))))\n",
    "        hess_inter = li@H@ri\n",
    "        H = hess_inter + (r*((s@(s.T)))) # BFGS Update\n",
    "        nabla = nabla_new[:] \n",
    "        x = x_new[:]\n",
    "        print(x)\n",
    "        it += 1\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-4.8  2.4]\n",
      "[ 4.8 -2.4]\n",
      "0.25\n",
      "[-1.57374114e-11  4.00000000e-01]\n",
      "[-0.06390533  0.49349112]\n",
      "1\n",
      "[-0.06390533  0.89349112]\n",
      "[0.05149931 0.09268051]\n",
      "1\n",
      "[-0.01240602  0.98617163]\n",
      "[0.01219207 0.01388197]\n",
      "1\n",
      "[-2.13951443e-04  1.00005361e+00]\n",
      "[ 2.12410723e-04 -4.40313684e-05]\n",
      "1\n",
      "[-1.54072053e-06  1.00000957e+00]\n",
      "[ 1.58870417e-06 -9.48722823e-06]\n",
      "1\n",
      "[4.79836375e-08 1.00000009e+00]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([4.79836375e-08, 1.00000009e+00])"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_opt = BFGS(f,[-1.2,1],100)\n",
    "x_opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bfgs_newton(f, x, iters):\n",
    "    \"\"\"\n",
    "    实现BFGS拟牛顿法\n",
    "    :param f: 原函数\n",
    "    :param x: 初始值\n",
    "    :param iters: 遍历的最大epoch\n",
    "    :return: 最终更新完毕的x值\n",
    "    \"\"\"\n",
    "    # 步长。设为1才能收敛，小于1不能收敛\n",
    "    learning_rate = 1\n",
    "    # 初始化B正定矩阵\n",
    "    B = np.eye(2)\n",
    "    x_len = x.shape[0]\n",
    "    # 一阶导g的第二范式的最小值（阈值）\n",
    "    epsilon = 1e-5\n",
    "    for i in range(1, iters):\n",
    "        g = grad(f, x)\n",
    "        if np.linalg.norm(g) < epsilon:\n",
    "            break\n",
    "        p = np.linalg.solve(B, g)\n",
    "        print(p)\n",
    "        # 更新x值\n",
    "        x_new = x - p*learning_rate\n",
    "        print(\"第\" + str(i) + \"次迭代后的结果为:\", x_new)\n",
    "        g_new = grad(f, x_new)\n",
    "        y = g_new - g\n",
    "        k = x_new - x\n",
    "        y_t = y.reshape([x_len, 1])\n",
    "        Bk = np.dot(B, k)\n",
    "        k_t_B = np.dot(k, B)\n",
    "        kBk = np.dot(np.dot(k, B), k)\n",
    "        # 更新B正定矩阵。完全按照公式来计算\n",
    "        B = B + y_t*y/np.dot(y, k) - Bk.reshape([x_len, 1]) * k_t_B / kBk\n",
    "        x = x_new\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[      0.         1403693.16074048]\n",
      "第1次迭代后的结果为: [ 1.00000000e+01 -1.40367316e+06]\n",
      "[ 935799.92023439 -311923.60062558]\n",
      "第2次迭代后的结果为: [ -935789.92023439 -1091749.5601149 ]\n",
      "[ -908836.82207169 -1044568.87169601]\n",
      "第3次迭代后的结果为: [-26953.0981627 -47180.6884189]\n",
      "[-26975.04462235 -47060.69516394]\n",
      "第4次迭代后的结果为: [  21.94645965 -119.99325495]\n",
      "[  21.40970748 -120.80932288]\n",
      "第5次迭代后的结果为: [0.53675217 0.81606793]\n",
      "[ 0.53609832 -0.18468436]\n",
      "第6次迭代后的结果为: [6.53854772e-04 1.00075229e+00]\n",
      "[0.00065329 0.0007513 ]\n",
      "第7次迭代后的结果为: [5.68951527e-07 1.00000099e+00]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([5.68951527e-07, 1.00000099e+00])"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bfgs_newton(f, np.array([10, 20]), 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dfp_newton(f, x, iters):\n",
    "    \"\"\"\n",
    "    实现DFP拟牛顿算法\n",
    "    :param f: 原函数\n",
    "    :param x: 初始值\n",
    "    :param iters: 遍历的最大epoch\n",
    "    :return: 最终更新完毕的x值\n",
    "    \"\"\"\n",
    "    # 步长\n",
    "    learning_rate = 1\n",
    "    # 初始化B正定矩阵\n",
    "    G = np.eye(2)\n",
    "    x_len = x.shape[0]\n",
    "    # 一阶导g的第二范式的最小值（阈值）\n",
    "    for i in range(1, iters):\n",
    "        g = grad(f, x)\n",
    "        if np.linalg.norm(g) < 1e-5:\n",
    "            break\n",
    "        # 更新x值\n",
    "        # print(G, g)\n",
    "        # H = \n",
    "        # print(H)\n",
    "        # learning_rate = line_search(f,x,H,g) \n",
    "        x_new = x - G@g* learning_rate\n",
    "        print(\"第\" + str(i) + \"次迭代后的结果为:\", x_new)\n",
    "        g_new = grad(f, x_new)\n",
    "        y = g_new - g\n",
    "        k = x_new - x\n",
    "        # 更新G正定矩阵\n",
    "        G = G + k.reshape([x_len, 1]) * k /(k@y)- (G@y).reshape([x_len, 1]) * (y@G)/ (y@G@y)\n",
    "        x = x_new\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "第1次迭代后的结果为: [1.00000000e+00 8.25721859e+04]\n",
      "第2次迭代后的结果为: [38110.13019904 52931.15030355]\n",
      "第3次迭代后的结果为: [11167.37988982 17450.01716326]\n",
      "第4次迭代后的结果为: [-14.94960005  33.18743898]\n",
      "第5次迭代后的结果为: [-0.42871266  1.34024197]\n",
      "第6次迭代后的结果为: [-1.31327403e-04  9.99827546e-01]\n",
      "第7次迭代后的结果为: [-8.65826032e-08  9.99999865e-01]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([-8.65826032e-08,  9.99999865e-01])"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dfp_newton(f, np.array([1, 2]), 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "10e-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "第1次迭代后的结果为: [1.00000000e+00 8.25721859e+04]\n",
      "第2次迭代后的结果为: [55047.09066136 64222.31370208]\n",
      "第3次迭代后的结果为: [1585.36647864 2775.41166488]\n",
      "第4次迭代后的结果为: [-1.25613672  8.15953999]\n",
      "第5次迭代后的结果为: [-0.03144786  1.0108145 ]\n",
      "第6次迭代后的结果为: [-3.91270926e-05  9.99955150e-01]\n",
      "第7次迭代后的结果为: [-3.42715061e-08  9.99999940e-01]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([-3.42715061e-08,  9.99999940e-01])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.15"
  },
  "vscode": {
   "interpreter": {
    "hash": "d99a3f7b344b3c3107482760db15f42178bfad658d282ab0a919b76809e13cb5"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
